Quantum Monte Carlo study of S = 1/2 weakly-anisotropic antiferromagnets on the 

square lattice 
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We study the finite-temperature behaviour of two-dimensional 5 = 1/2 Heisenberg antiferromag- 
nets with very weak easy-axis and easy-plane exchange anisotropies. By means of quantum Monte 
Carlo simulations, based on the continuous-time loop and worm algorithm, we obtain a rich set of 
data that allows us to draw conclusions about both the existence and the type of finite-temperature 
transition expected in the considered models. We observe that the essential features of the Ising 
universality class, as well as those of the Berezinskii-Kosterlitz-Thouless (BKT) one, are preserved 
even for anisotropies as small as 10""^ times the exchange integral; such outcome, being referred 
to the most quantum case S = 1/2, rules out the possibility for quantum fiuctuations to destroy 
the long or quasi-long range order, whose onset is responsible for the Ising or BKT transition, no 
matter how small the anisotropy. Besides this general issue, we use our results to extract, out of the 
isotropic component, the features which are peculiar to weakly anisotropic models, with particular 
attention for the temperature region immediately above the transition. By this analysis we aim to 
give a handy tool for understanding the experimental data relative to those real compounds whose 
anisotropies are too weak for a qualitative description to accomplish the goal of singling out the 
genuinely two-dimensional critical behaviour. 

PACS numbers: 75.10.Jm, 75.40.-s, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 



In the last few decades the Heisenberg antiferromag- 
net (HAFM) on the square lattice has been thoroughly 
studied by means of seiieml theoretical, numerical and ex- 
perimental techniquesBBu. Such research hands us a pic- 
ture where three classes of substantially different models 
appear, the isotropic, the easy-axis, and the easy-plane 
antiferromagnets. The reference models for the two lat- 
ter classes are the 2D-Ising and 2D-XY antiferromagnets, 
defining their respective universality class and known to 
display a finite-temperature phase transition. It is also 
known that quantum effects, despite being enhanced by 
the reduced dimensionality, do not substantially change 
the essential features of the above models, not even in 
the extreme S = 1/2 case, though evidently renormal- 
izing their thermodynamic properties. This results frogi 
renormalization groupH'El and semiclassical approachesBQ, 
as welL jas_&xiin quantum Monte Carlo (QMC) simula- 
tionfflBEiy. 

In the above picture, however, there still is a grey area, 
where our knowledge is not detailed enough to allow a 
precise reading of the experimental data: this is the area 
of very weak anisotropies and strong quantum effects, 
which is of particular interest as most of the real lay- 
ered compounds whose magnetic behaviour is properly 
modeled by the 5=1/2 HAFM on the square lattice 



are characterized by anisotropies as small as 10^"^ times 
the exchange integral. These compounds exhibit a phase 
transition at a finite temperature Tn which is often too 
large for the interlayer coupling to be the unique player, 
while the idea of a two-dimensional anisotropic critioaJi 
as the trigger of the transition appears well soundtJl 
such idea is corroborated by the measured values of some 
critical exponentsll3. The experimental observation tells 
us that three-dimensional long-range order is present be- 
low Tn and that well above Tn no trace of anisotropic 
behaviour is left; it is slightly above Tn that one hence ex- 
pects evidences of genuine 2D anisotropic behaviour to be 
detectable. In order to let these experimental evidences 
surface out of the sea of the isotropic thermodynamics, 
precise numerical data for the 5 = 1/2 nearly isotropic 
HAFM are needed: it is the purpose of this work to fulfill 
such need. In particular, we address the problem of the 
existence and characterization of the phase transition, 
making use of continuous-time QMC cluster algorithms, 
which are well suited in the neighbourhood of a critical 
point, since they do not suffer from critical slowing-down. 
The general resulting picture is that a phase transition is 
induced by an arbitrary amount of anisotropy, and that 
several distinctive features of the expected universality 
class can be traced out. 

The structure of the paper is as follows: in Section || 
the general properties of the anisotropic models are de- 
scribed, extracting the most significant difference with 
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respect to the isotropic model. In Section III the QMC 
methods are presented, the thermodynamic quantities 
under investigation are defined, and the finite-size scal- 
ing (FSS) theory used in our analysis is briefly recalled. 
The results for nearly isotropic models, as from both FSS 
analysis and thermodynamic behaviour, are presented 
and discussed in Sections IV and ^ for the easy-axis and 
easy-plane case, respectively. In Section VI the critical- 
temperature vs anisotropy phase-diagram is discussed. 
Eventually, conclusions are drawn in Section VII. 



II. 



THE TWO-DIMENSIONAL XXZ MODEL ON 
THE SQUARE LATTICE 



The XXZ model is defined by the Hamiltonian 

J • 

+ (1-Aa) , (1) 



where i = (ii,?2) runs over the sites of an L x L square 
lattice, d connects each site to its four nearest neigh- 
bours, J > is the antiferromagnetic exchange inte- 
gral, and are the easy-axis (EA) and easy-plane 
(EP) anisotropy parameters, respectively, hereafter get- 
ting positive values smaller than unity. As J sets the 
energy scale, the dimensionless temperature t = k^T / J 
will be used in the following. 

The spin operators Sf {a = x, y, z) obey the sm(2) 
commutation relations [5*1*, S*^] = ie'^^'^SijSj and are 
such that ISp =5'(5'-hl). 

Besides the isotropic model A^ — A^ — 0, Eq. (|l|) 
defines the EA (Aa = 0, < A^ < 1) and EP (A^ = 0, 
< Aa < 1) magnets, whose respective reference models 
are the Ising model (A^ — 1) and the XY (also known 
as XXO) one (Aa = 1). 

What is known about these models can be summarized 
as follows: 

i) the isotropic model has no finite-temperature tran- 
sition; its ground state is ordered for any S and a critical 
region of divergent correlations is clearly observed at very 
low temperatures; 

ii) the EA models exhibit an Ising-like transition at a 
critical temperature t-^ which is an increasing function of 
both A^ and S*; for S* > 1, is finite for all anisotropics; 

iii) the EP models exhibit a transition of the 
Berezinskii-Kosterlitz-Thouless (BKT) type, at a critical 
temperature ^b^t which is an increasing function of both 
A^ and S] for 5 > 1, ^b^t finite for all anisotropics. 

Some of the above statements are rigorously proved, 
others come from the combination of theoretical, numer- 
ical and experimental results. 

More precisely, for what concerns ground state prop- 
erties, an antiferromjametically ordered ground state has 
been exactly provedtatZl for 5 = 1/2 only in the case of 
strong anisotropy, Aa > 0.88 and A^ > 0.32. For weaker 



anisotropics clear numerical evidence exists of an ordered 
ground state down to the isotropic limit, .as coming from 
QMC and exact diagonalization studieala. As for the 
finite-temperature jncoperties, according to renormaliza- 
tion group analysia3't3, the critical temperature vanishes 
in the isotropic limit as ^i^bkt ^ [ln(const./A^^A)]^^- In 
the EA case the existence of the transition for arbitrary 
A^ has been rigorously proved only in the classical limit 
5 — > oo, while in the quantum case the proof is reatricted 
to a spin-dependent finite amount of anisotropjCj . In 
the EP case the situation is even less clear, since also in 
the classical limit a rigorous proofcil of the existence of 
a transition is limited to the case Aa ^ 0.553, while no 
rigorous proof exists in the quantum case. 

In the S = 1/2 case, evidences of a phase transi- 
tionJbr anisotropics as small as A^ ^p-O.Ol in the EA 
caseaci and Aa — 0.02 in the EP caseEHl were suggested 
by previous QMC approaches; however, these computa- 
tions employed local algorithms which cannot easily ac- 
cess the critical region and a rigorous FSS analysis could 
not be performed. The situation is still unclear also be- 
cause r«jGent works based on real-space renormalization- 
groupE3'E3 predict the existence of a critical value of the 
anisotropy in the EA case (A^'^^ f« 0.2) below which the 
transition would be destroyed by quantum fiuctuations. 
In this work we hence consider 5* = 1/2 and four nearly 
isotropic systems, two EA, A^ 
and two EP ones, Aa = 0.02 and Aa = 0.001 



0.01 and A^ = 0.001, 



III. QUANTUM MONTE CARLO, 
OBSERVABLES AND FINITE-SIZE EFFECTS 

A. Quantum Monte Carlo method: 
continuous-time algorithms 

As usually done in the existing literature on QMC, in 
this section and in appendices ^ and ^ we employ the 
notation 



J 



XY 



= J(l - Aa) 



(2) 



and (3 = l/keT. 

The QMC method for the 5=1/2 XXZ model is 
based on the Trotter-Suzuki decomposition of the par- 
tition function, whfch can be approximated by the fol- 
lowing expressiorcJ: 



Z{f3) = Tr e' 



En 



(Ar) 



(3) 



where Wp represents the amplitude of propagation of 
a pair of nearest-neighbour spins from a configuration 
\ui^Uj) to \cf'i,<y'j) in the (imaginary-) time step At = 
(i/M, M being the Trotter number and \{<Ji}) (cTj = 
±1/2) the basis set diagonalizing the 5*1 operator. The 
two bond configurations define a space-time plaquette 
configuration p = {ci, (Tj; cr^, (7^}, so that Wp can be seen 
also as the weight of a given plaquette configuration p. 



(1) 





(2) 



FIG. 1: Plaquette configurations with non-vanishing weights 
in the 5 = 1/2 XXZ model. The vertical axis is the imaginary- 
time direction. 



The index n runs over all plaquettes on the space-time 
lattice, and the index S runs over all configurations of 
the system. At each time step plaquettes are defined 
on different groups of bonds (ij) , so that all bonds in- 
volved in the propagation at the same time step do not 
share any spin; moreover, each plaquette shares its corner 
spins with two plaquettes on the previous and two on the 
subsequent time step. The expression (||) becomes exact 
in the limit M ^ oo. 

In the case of the XXZ model, only the 6 plaquette 
configurations shown in Fig. ^ have non-zero weights 
whose expansion to first non trivial order in At are"^^" 
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(4) 



Plaquettes 1, 2, 3, and 4 propagate the state of the spin 
pair unchanged, while plaquettes 5 and 6 introduce an 
exchange of state for the interacting spins, hereafter de- 
noted as kink. 

Any standard Monte Carlo (MC) method can be 
applied for numerically evaluating the partition func- 
tion (H) ; in principle one could do single local Metropolis 
moves which update clusters of plaquettes sharing cor- 
ner spins in_prder to move to configurations with non- 
zero weighted. This approach is however inconvenient, 
being affected by the well-known critical slowing-down. 
Moreover, in order to control the Trotter-approximation 
error, it is necessary to repeat the simulation for in- 
creasing Trotter number. In recent years efficient clus- 
ter algorithms were proposed to overcome the problem 
of critical slowing down in QMC simulations. In the 
context of quantum spin systems, the loop-cluster algo- 
rithm and the worm algorithm have revealed very effec- 
tive; moreover, they can be also formulated in continuous 



G: 



FIG. 2; Breakup decompositions of a single plaquette in the 
5=1/2 XXZ model; thick lines join grouped spins. 



imaginary time, thus removing completely the Trotter- 
approximation error. 



1. Loop algorithm 

One can completely eliminate the critical slowing-down 
by introducing the so-called loop-cluster algorithrnEBj. 



which is tie quantum analog of the Swendsen-WangE 
and Wolffn^l cluster algorithms introduced for classical 
spin systems. Within the multi-cluster approachlij (ana- 
logue to the Swendsen-Wang scheme) the loop algorithm 
consists of probabilistically assigning to each plaquette 
a break-up decomposition (or graph) G, i.e., a way of 
grouping its spins in subgroups, so that the grouped spins 
can be flipped all at once bringing the plaquette into a 
configuration with non- vanishing weight; in the case of 
the XXZ model the above condition allows only group- 
ing of spins into pairs (non-freezing breakups) or all to- 
gether (freezing breakup), as shown in Fig. g. Assigning 
a breakup G to each plaquette, one univocally defines a 
breakup decomposition Q = {Gn} of the whole configu- 
ration S = {(Ti^r} into loops. Different breakup decom- 
positions G are assigned to each plaquette configuration 
p according to weights w(p, G) obeying the general sum 
rule Wp = w(p, G). Since the XXZ model is symmet- 
ric under time reversal, it is straightforward to require 
that plaquettes connected by time-reversal must have the 
same weights, i.e., w(l,G) = w(2,G), w(3,G) = w(A,G) 
and w(5,G) — w(6,G). These weights must be positive 
and obey the detailed balance condition. The presence 
of frozen plaquettes leads to the generation of very big 
loops, whose flip is expected to be less and less effec- 
tive in generating successive uncorrelated configurations; 
therefore, a general strategy to be followed is that of min- 
imizing the probability of frozen plaquettes to appear. If 
one sets ^(p, (g)) = for each p, a solution with positive 
breakup weights is obtained only in the case J^^ > , 
i.e., in the EP and isotropic case. When > J'^^, 
i.e., in the EA case, one must allow for freezing of at 
least a single plaquette configuration to achieve positive- 
ness. In particular, in the antiferromagnetic case, it is 
w(3,(X)) = w(4,(g)) that must be different from zero in 
order to have positive weights. By further imposing in- 
dependence of the breakup decomposition from the pla- 
quette configuration, i.e.. 



w(l,\\)^w(3., 
w(3,= 



w(l, x) 
w(5,=) 



w(5, x) 



(5) 
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one gets that, within the same loop decomposition Q, all 
of the configurations S' obtained from the starting one S 
by flipping some of the loops, are reached with the same 
probability, i.e., each loop is flipped independently of all 
the other ones. 

The complete set of non-zero breakup weights used in 
our calculation thus readsE3 (to first order in At): 



• Easy-plane and isotropic case {Jx > Jz)- 



1 —At 



J 



XY 



x) 

w;(3,=) 
• Easy-axis case {Jx < Jz)- 
^1,11) 



J^^ - J' 



■At 



-At 



(6) 



w(3,«)) 
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XY 



-At 



(7) 



The fact that freezing is required in the EA case can 
be physically interpreted as a reflection of the enforced 
antifcrromagnctic correlations of the z-component of the 
spins: as a matter of fact, freezing of plaquettes of type 
3 and 4 implies preservation of nearest-neighbour anti- 
ferromagnetic correlation along the MC flow, i.e., anti- 
ferromagnetic configurations are more "resistant" to MC 
fluctuations, with respect to ferromagnetic onps. 

In the limit of continuous imaginary timetj. At — ^ 0, 
plaquettes with no kink, i.e., of types 1 (2) and 3 (4), ac- 
quire unitary weight, while plaquettes with a kink, i.e., 
of type 5 (6), get an inflnitesimal weight, still keeping a 
finite weight per unit time Wp = limAr^o u^p/At; there- 
fore kink-bearing plaquettes must be regarded as Pois- 
sonian events in the imaginary-time evolution of each 
pair of interacting spins. At the same time the breakup 
decomposition creating no kink in imaginary time evo- 
lution, acquires a unitary weight, while all the other 
breakups acquire infinitesimal weight, still keeping a fi- 
nite weight per unit time w(p, G) . In the case of plaquette 
5 (6) , since the breakup weights have to be normalized to 
the plaquette weights, they become finite probabilities 



M5, 



=) = ^l + jS-) M5,x) = l-p(5,=) (8) 



in the EP case, and ^(5,=) = 1, p(5, x) = in the EA 
case. 

The algorithm then proceeds as follows: 

• distribute breakups x , = and (8) on the continu- 
ous segments (= infinite sequences of infinitesimal 



plaquettes 1(2) and 3(4)) along the imaginary-time 
evolution of each pair of interacting spins, accord- 
ing to the Poisson distribution having as parameter 
/3w(p, G); for each kink in the propagation, choose 
a breakup with probabilities p(5, G). 

• reconstruct the loops defined by the decomposition 
of each infinitesimal plaquette; 

• decide whether to fiip each loop independently with 
probability one half. 

The above procedure (multi-cluster update) represents 
a single MC step in our code. We have generally per- 
formed 10* MC steps for thermalization for each value 
of the temperature, and 1 1.5 • 10^ MC steps for eval- 
uation of thermodynamic observables. The algorithm is 
very efficient in both the EA and the EP case, with auto- 
correlation times which always remain around unity for 
all the lattice sizes L we considered, i.e., L —16, 32, 64, 
128, 200. The autocorrelation time Tc has been estimated 
by the blocking techniqu£3 as: 



(9) 



where cr^ is the variance of the time-series {xi} 
(i = 1, A'stops = n-bA'b) produced for the variable x, 
while tr^ is the variance of the block variable Xj = 

^b"' ES1i)*jv,+i ^^ ij = l,-:«b), with » T for 
the estimate to be sensible. 

The introduction of freezing breakups has the effect of 
making a loop branch into subloops; as the flip of large 
branched clusters is less effective in generating succes- 
sive uncorrelated conflgurations, the inclusion of freezing 
is generally thought to lower the efficiency of the loop al- 
gorithrrEj, although no direct evidence of such conclusion 
exists. Nevertheless, for the EA anisotropics we consider, 
no significant loss of efficiency (i.e., increase in correlation 
time) is observed. Moreover, in the Ising limit, the above 
algorithm reproduces the Swendsen-Wang algorithm for 
the Ising model, which is known not to suffer from slow- 
ing down; therefore efficiency is possibly lost only for 
intermediate values of the anisotropy, which is not the 
case we are interested in. rin 

We have implemented improved estimatorsEj^t2l for all 
the quantities of interest. A separate, more careful analy- 
sis is needed in the case of off-diagonal observables, whose 
most general bihnear example may be (S^ {t)SJ {t')) . 
In ab stjvne e of freezing, the improved estimator simply 
readsy^: (5+(t)^7 (t'));^^^^ = 1, ff {i,r) and {j,r') 
belong to the same loop, = otherwise. This result 
can be understood by observing that the quantity we 
are considering may be seen as the ratio between the 
partition function of the system with density opera- 
tor Sf{T)Sj{T')eyi\){—(3'H) and the partition function 
of the original model considered, with density operator 
exp(— /37i). The final outcome for the estimator given 
above thus follows from the observation that in absence 
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of freezing, i.e., in the isotropic or EP class, according to 
Eqs. (|^) the weight assigned to a loop decomposition is 
actually independent of the configuration onto which it is 
defined, since each kind of breakup has the same weight 
on every plaquette configuration where its weight is non- 
zero; therefore the same loop decomposition is reached 
with equal probability in the original and in the modi- 
fied system, with the additional constraint, in the case of 
the modified system, that the two points {i, r) and {j, t') 
must belong to the same loop. 

When freezing is present, the above argument breaks 
down, since only one plaquette configuration admits 
freezing and the constraint of having (i, r) and (j, t') on 
the same loop is no longer sufficient to have a non-zero 
contribution to the estimator. In principle it is possible 
to define the estimator for off-diagonal observables even 
in case of freezing; however, not only its implementation 
is highly non trivial from the point of view of program- 
ming, but its evaluation would also consume a consid- 
erable amount of computational time. As a matter of 
fact, for each pair of points (i,T), (j,t') belonging to 
the same branch of a loop, one should look for freezing 
breakups by checking all the loop segments between the 
two points, an operation that, being the loops structure 
involved, becomes very costly at low temperatures. We 
have therefore renounced to implement such estimators 
in case of freezing. To have a complete picture of the ther- 
modynamics of the system in the EA case we have then 
resorted to a different (and generally less efficient) QMC 
scheme, within which the calculation of the off-diagonal 
observables in the EA case is relatively straightforward, 
i.e., to the so-called worm algorithm. 



2. Worm algorithm 

The worm algorithm represents an alternative way to 
overcome the problem of critical slowing down in QMC 
simulations. The original idea of the algorithm can be 
found in Ref. but we here formulate the algorithm in 
a different way, so that it appears as a direct generaliza- 
tion of the loop algorithm; our formulation is more di- 
rectly related with the so-called "operator-loop update" 
introduced in the framework of the stochastic series ex- 
pansionEa. 

The worm algorithm starts by choosing a point at ran- 
dom in space-time, inserting two discontinuities in the lo- 
cal imaginary-time evolution, and then keeping one fixed 
(the "tail" of the worm) while letting the other (the 
"head" of the worm) freely travel through the lattice. 
The single-worm update ends when the head happens 
to "eat" the tail (the worm closes), so that the isolated 
discontinuities disappear and the system is led to a new 
configuration having non-zero weight. All the segments 
of imaginary-time evolution touched by the worm's head 
have to be fiipped, i.e., the worm's head performs a real- 
time update of the system. Its motion conventionally 
goes forward (backward) in imaginary time while updat- 



ing segments with up (down) spins, and it is ruled by 
detailed balance condition, to be locally satisfied on each 
(infinitesimal) plaquette it touches. In the context of the 
worm algorithm one has to abandon the concept of pla- 
quette breakup, and the decision about how a plaquette 
is passed through by the worm must be taken from time 
to time. One of the distinctive feature of the algorithm 
is that the worm can pass through the same plaquette 
many times, each time in a different way, depending on 
the local configuration left at its previous passage. 

General detailed balance conditions for the single pla- 
quette update when the worm's head passes through it, 
hipping two spins, read 



(10) 



where we have already introduced the time- and space- 
reversal symmetries, so that here "1" means 1 or 2, "3" 
means 3 or 4, and " 5" means 5 or 6, depending on the way 
the worm's head travels through the plaquette. More- 
over, the transition probabilities must satisfy the sum 
rules 



wi p{l- 


-3) 


= W3 p{3- 
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Wi p{l- 


-5) 


= W5 p{5- 




W3 p{3- 


-5) 


= W5 p{5~ 
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p{l- 


-3)- 


^P{1- 


-.5) = 
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p{3- 






-.5) = 


1 


p{5- 


-1)- 


^p{5- 


-.3) = 


1 



(11) 



Moving to the continuous-time limit, we express the tran- 
sition probabilities from a plaquette with a non- vanishing 
weight to another plaquette, as 



p(1^3) 
P(3^1) 
p(1^5) 
p(3->5) 



1 + 7r(1^3) At 
1 + 7r(3^1) Ar 
7r(l-^5) Ar 
7r(3^5) Ar , 



(12) 



where 7r(l-^5) and 7r(3^5) have the meaning of (pos- 
itive) transition probabilities per unit imaginary time, 
while 7r(1^3) and 7r(3^1) are (negative) corrections to 
the transition probabilities among plaquettes taking non- 
vanishing weights; these corrections arise from the pois- 
sonian occurrence of kinks in imaginary-time evolution. 
With the above parameterization of the transition prob- 
abilities, the first two of the sum rules (O) take the form 



7r(l- 
7r(3- 



.3) 
>1) 



7r(l- 
7r(3- 



>5) ^ 
>5) = 



(13) 



while the third remains unchanged, as p(5— >1) and 
p(5— >3) keep their meaning of dimensionless probabilities 
for the different ways the worm's head can pass through 
a kink in imaginary-time evolution. The set of detailed 
balance equations in the continuous-time limit takes the 
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form 



7r(3- 


-1) = 
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7r(l- 


-.5) = 


2 






2 


7r(3- 


-.5) = 



p(5^1) 

p(5--3) . (14) 



Together with the sum rules, they give as unique solution 
the following set of transition probabilities: 



7r(1^5) 
7r(3^5) 
p(5^3) 
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(15) 



It is immediate to see that this solution is equivalent to 
the set of breakup weights (^), which means that, at this 
level, the worm algorithm is nothing but the Wolff-type 
(single-cluster) version of the loop algorithm. However, 
as seen in the loop algorithm, in the EA case transi- 
tion probabilities are not always positive, and some other 
transition mechanism must be invoked to overcome this 
problem. As seen before, the remedy in the case of the 
loop algorithm was to allow for branching of the loops; 
if one hence allows for branching also in the worm algo- 
rithm, the single-cluster version of the loop algorithm for 
the EA case is obtained. 

A different strategy can be adopted in the case of the 
worm algorithm by introducing a new type of motion, 
named bouncing^ where the worm's head, when attempt- 
ing to update a plaquette, is bounced off and hence forced 
to locally trace back its route. From the physical point of 
view, the existence of a bounce mechanism protects some 
plaquettes from being updated, possibly those plaquettes 
containing local spin configurations which give a relevant 
contribution to the thermodynamics of the system; such 
local configurations are preserved along the MC flow and 
the effect of bouncing is seen to be very similar to that of 
freezing in the context of the loop algorithm. In the case 
of the EA antiferromagnet, the most relevant local con- 
figurations are those containing antiferromagnetic corre- 
lations of the z-components, i.e., in terms of plaquettes, 
p=3,4. Therefore, we allow for bounce on these plaque- 
tte configurations, introducing a finite bounce probability 
p(3, h) = 7r(3, 5) At which has to be accounted for in the 
sum rule 



7r(3^1) + 7r(3^5) + 7r(3, 6) = 



(16) 



The detailed balance condition for the bounce probabi lity 
is trivial, reading p( 3, 5) u;(3) = p(3, &) w(3). Eqs. ( |l5| ) 
and (|l^) form an underdimensioned set and 7r(3, b) can be 
hence chosen arbitrarily, with the only constraint of pos- 
itive transition probabilities. As in the case of freezing. 



it is highly convenient to minimize the bounce probabil- 
ity: when the worm's head bounces, part of its update 
operations are lost as it locally traces back its way, so 
that the efficiency in updating the configuration, keep- 
ing the number of elementary update operations fixed, is 
lowered. The following solution for the transition proba- 
bilities, minimizing the bounce probability, is found 



^(l->5) 
7r(3^5) 

^(3,5) 
p(5^1) 




2 

1 . 



tXY\ 



(17) 



The worm algorithm with the bounce process is a pure- 
quantum cluster algorithm: in the Ising model, which is 
a substantially classical statistical model, the algorithm 
loses its cluster nature, since only bounce processes sur- 
vive, thus confining the worm on a single site. In this case 
the worm algorithm reproduces the single-flip Metropolis 
algorithm, as shown in Appendix 

As in the case of the loop algorithm, each of our sim- 
ulations consists of lO'* MC steps for thermalization and 
of 1 -f- 1.5 * 10^ MC steps for evaluation of thermody- 
namic observables. During thermalization, the number 
of worms to be produced at each step is adjusted so that 
the total length of the worms in the imaginary-time di- 
rection roughly equals the size of the (D-l-l)-dimensional 
lattice, * (3; this number is then kept fixed during the 
measurement phase. In this way, autocorrelation times 
of the order of unity are achieved for all values of the 
EA anisotropy considered. At variance with the loop 
algorithm, the efficiency is here expected to drastically 
decrease as the anisotropy increases, given that, as the 
model moves toward the Ising limit, the cluster algorithm 
transforms into the local Metropolis algorithm; however, 
the case of strong anisotropy is not of our interest here. 

As mentioned in the previous section, the estimator for 
bilinear off-diagonal quantities like {t)SJ [t') can be 
thought of as a partition function for a modified model, 
in which two spin discontinuities are inserted in the sys- 
tem configuration at the points (i, r), (j, r'). Now it be- 
comes clear that configurations giving a non-zero contri- 
bution to such a partition function are generated during 
the worm update whenever the discontinuities associated 
to the head and tail of the worm coincide with the above 
points, both in the EA and in the EP case. Therefore 
the off-diagonal observables are_measured on-the-fly dur- 
ing the motion of worm's headEj, and each worm update 
produces a statistics for the estimators which grows lin- 
early with the length of the worm. As in the case of 
the loop algorithm, this kind of estimators give to the 
worm geometry a physical meaning: the further the worm 
head travels away from its tail, the larger will be the off- 
diagonal correlations. On the other hand, improved esti- 
mators are not defined for diagonal quantities in the EA 
case; to this respect, worm and loop algorithms are seen 
to be exactly complementary. 
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As final remarks, we mention that the worm algorithm 
retains its full efficiency also in presence of a uniform 
magnetic field applied to the spins, while the loop algo- 
rithm is known to exponentially slow down as-the field is 
increased and/or the temperature is lowerecO. Details 
of the worm algorithm in a finite field are not relevant 
to the present work and will be given elsewhere. Inde- 
pendently of us, Syljuasen and SandvikEZI have recently 
developed a very similar (directed loop) algorithm, in the 
framework of both stochastic series expansion and path- 
integral Monte Carlo. 



B. Thermodynamic quantities 

We briefly report here the definition of the relevant 
thermodynamic quantities measured in our QMC study, 
together with their respective estimators. The MC aver- 
age of the estimator will be hereafter denoted as {...)^^. 

The internal energy (7i) is estimated as the MC aver- 
age of 



i.d 



dr 0i,d(T) = E 



(18) 



where 4>i^d{T) takes the value — if at imaginary time 
T there is an infinitesimal plaquettc configuration of type 
1 (2), J^/4ifoftypc3 (4), -6{t) if of type 5 or (6). This 
corresponds to the continuous- time limit of the energy 
estimator as defined in Ref. [27[ 

The specific heat c = P^HH^) - {'H)^)/L'^ is estimated 
from energy fluctuations as 



1 

l2 



((/3'£^'-^kinks)„^-/3^(£;)L) 



(19) 



where TVkinks is the number of kinks present in each gen- 
erated conflguration. The variance of the speciflc heat 
has been estimated via binning analysis of the time se- 
ries related to the energy estimator and the kink number. 

The staggered magnetization = (^1)*('S'|) is esti- 
mated as the MC average of 

-l5](~l)V|EEm, . (20) 



The spin-spin correlation function is 

^""(^) drdT'{Sf{T)Sf+,{T')) /(r,r'), (21) 

where f{T,T') ~ (3 5{t — t') in the equal-time correlator 
(ET) and f{T,T') = 1 in the time-averaged (TA) one. 
In both cases the numerical calculation of the correlation 
function takes advantage of the existence of the improved 
estimator deflned in the previous subsection. 
The generalized susceptibility is 



(22) 



the time-averaged susceptibility corresponds to the ther- 
modynamic deflnition (second derivative of the free en- 
ergy) while the equal-time one corresponds to /3 * S{q), 
where S{q) is the static structure factor as measured, 
e.g., in neutron scattering experiments. From the general 
deflnition above follow those of the uniform susceptibility 



and of the staggered one 



(23) 



(24) 



Susceptibilities and correlation functions have been mea- 
sured both along the z-axis (C^^, x^^) ^'^d in the xy- 
plane (C^"= = C^^, x=""= = x"^)- in the EA case, the latter 
have been evaluated by means of the worm algorithm. In 
what follows, we will show and comment our data rela- 
tive to the uniform TA susceptibility and to the staggered 
ET susceptibility, being such quantities the more relevant 
ones from the experimental point of view. 

The correlation length is defined via the long- 
distance exponential decay of the staggered correlation 
function, (-l)''C""(r) - cxp(-r/f"") (r oo). A di- 
rect estimate of the correlation length may be hence 
found by fitting the long-distance behaviour of C""(r) 
with a model-dependent function, as discussed in the fol- 
lowing sections. Such procedure, however, is strongly de- 
pendent on the quality and stability of the fit, and does 
not always lead to a univocally defined result in case of 
a finite-size system in presence of a phase transition, i.e., 
of a diverging correlation length. An alternative strat- 
egy, which we have also used, is offered by the so-called 
second moment definitionE3 




X""(7r, tt) 



'{n + 2'K/L,n) 



1 



(25) 



which can be directly extracted by the simulation data, 
supplemented by a binning analysis of susceptibility time 
series in order to estimate the variance. 

Another relevant observable, in the EP case, is the 
helicity modulus T, which is a measure of the response of 
the system to the application of a twist $ in the boundary 
condition along a given direction: 



T 



(26) 



A=0 



where (j) = ^/L. In Appendix ^ we show that, starting 
from the above definition as explicitly written in terms 
of spin operators, the estimator of the helicity modulus 
of the 5* ^ 1/2 XXZ EP model reads 



t. 



IWI 



(27) 



where W — {Wi,W2), W^i(2) being the total winding 
number of spin paths (paths traced by a fixed spin con- 
figuration, up or down) in the 1(2) lattice direction. Re- 
markably, this estimator is directly related with that of 
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the superfluid density of bosonic systemgSa. An efficient 
improved version of the estimator, (|2^) has been intro- 
duced by Harada and Kawashimalla in the context of the 
loop algorithm, and is the one employed in this work. 



C. Finite-size scaling 

A FSS analysis^ can give strong indications on the ex- 
istence of a phase transition at some temperature tc , pos- 
sibly leading to a full characterization of its universality 
class. The simplest evidence that a transition occurs is 
found when, for increasing lattice size, the order param- 
eter scales to a finite value below a certain temperature, 
indicating that a non-zero order parameter develops in 
the thermodynamic limit. 

In the case of second-order phase transitions, the 
AnsatzcJ for the scaling behaviour of a generic finite-size 
thermodynamic quantity in the neighbourhood of 

the critical point reads 



AUt) ^ LPl^ Fa{L^I'' {t-t. 



(28) 



the 



critical 

p 



exponent of A 



where p 

A(t-^tc) ^ \t ~ tc| V is the exponent for the corre- 
lation length, while Fa is the universal scaling function. 
At the critical point Eq. (H) imphes A^it^) ~ LpI" . 
In the case of ^ this means linear scaling at critical- 
ity, without any assumption on the universality class; 
therefore, looking for the temperature at which a prop- 
erly defineccJ ^L(t) scales linearly with the system size 
gives an unbiased estimate of the critical temperature. 
Eq. ( Pq ) implies that the scaling plot of L~p/'^ vs y = 
{t — t^L^/'^ ^ with a proper estimate of <c, shows the data 
for different lattice sizes to collapse onto the universal 
curve Fa{]j)- 

In the case of a BKT transition, in which no order 
parameter is given, the presence of topological order at 
finite temperature is shown when the helicity modulus 
scales to a finite value below a certain temperature. The 
use of the scaling Ansatz to locate the critical tempera- 
ture can be generalized to the case of a BKT transition, 
though most of the critical exponents are not defined. 
However the Kosterlitz-Thouless theory predicts -q — \/A 
at the critical point, so that a scaling behaviour of the 
susceptibility as L'^~^ = Vl^ is a good signature of the 
critical temperature.pCyloreover, Kosterlitz's renormaliza- 
tion group equationsc3 provide tbe. critical scaling law for 
the helicity modulus in the formCJ 



1 



21og(L/Lo) 



(29) 



where Lq is a constant. This relation has been widely 
used to locate the BKT criticaLteinperature of the clas- 
sical 2D planar-rotator modelfja and of the S — 1/2 
quantum XY modelt2l. 

We end this section with a general remark. It is ob- 
served that the smaller the anisotropy, the bigger the 



lattice sizes required to enter the asymptotic scaling 
regime, where FSS holds. This is essentially due to the 
fact that the critical region is shifted to lower tempera- 
ture: the correlation length of the isotropic model, acting 
as a lower bound for that of the nearly isotropic ones, 
increases exponentially upon lowering the temperature, 
and therefore, keeping the lattice size fixed, the ratio 
L/^L, that drives the onset of asymptotic scaling near 
the transition, gets smaller. 



IV. EASY-AXIS MODEL AND ISING 
TRANSITION 



The values of the anisotropy here considered are = 
0.01 (also used in Ref. |) and = 0.001. They are 
comparable with the characteristic anisotropies of real 
compounds; yet, for such small anisotropy there is,, 
universal consensus onjihe existence of a transitionESi 
From previous workaJu'EI the transition is expected in 
the temperature range 0.2 < t < 0.3 in both systems. At 
higher temperature the behaviour gets closer to that of 
the isotropic model, which has been eaten s ivply investi- 
gated by means of QMC in recent yearaSJea'cZlH; we have 
extended our analysis up to t ~ 0.8 in order to identify 
those deviations from the isotropic behaviour that can 
be experimentally detected above the critical region. 

In our approach, evidences for the existence of an Ising- 
like transition follow from a detailed FSS analysis of the 
data; subsequently, we analyze the temperature depen- 
dence of some relevant thermodynamic quantities, em- 
phasizing the signatures of the EA nature. 



A. Finite-size scaling analysis 

Our analysis proceeds in three steps: we give evidence 
for a transition to occur, then the transition tempera- 
ture is located, and eventually the Ising critic al sca ling 
is tested. After the discussion made in Section HI C the 



FSS analysis for A^ = 0.001 is expected to be more deli- 
cate than for A^ ~ 0.01. Indeed, for the lattice sizes used 
(L < 128) some quantities show to have well entered the 
asymptotic scaling regime, while others do not. Anyway, 
clear (though not complete) evidence for the Ising uni- 
versality class is given also for A^ = 0.001; larger lattices 
would be required to reach a full characterization. 

Let us first consider the order parameter, i.e., the stag- 
gered magnetization given in Eq. (|2^) . In Fig. || ill for 
A^ = 0.001 is seen to scale to a finite value ii t ^ 0.22, so 
that the magnetization in the thermodynamic limit be- 
comes finite; the same behaviour is a fortiori observed 
in the case A^ = 0.01. We then invoke the scaling 
Ansatz (|2^) for the longitudinal correlation length 
The scalin g plot of as specifically defined in Sec- 
is shown in Fig. 
=0.001) = 0.2225(15) 
=0.01) = 0.2815(25). 



tion 



IV E below. 



and gives ii(A^= 
ysis yields ii(A^: 



for A,, 



0.001 



A similar anal- 
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FIG. 3: Scaling of the staggered magnetization in the EA 
model with = 0.001, for different t. 
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FIG. 5: Scaling of the Binder's fourth cumulant in the EA 
model with — 0.01 for different t. The solid line indicates 



the universal critical value ' (see text) 
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FIG. 4: Scaling of the longitudinal correlation length in 
the EA model with A^ = 0.001, for different t. 



Hitherto, no assumptions were made about the uni- 
versality class. In order to identi% it, we consider the so 
called Binder's fourth cumulantll3, shown in Fig. ^ and 
defined by 



U4 ^ 1 



\ a ' MC 



(30) 
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lattice 



00(1) at in the 2D Ising model on the square 
J, and increases (decreases) with L, below (above) 
tj. For = 0.01, we verify such behaviour and ob- 
tain tj(A^=0.01) = 0.280(3), consistently with the above 
unbiased estimate from the scaling of f^^. The scal- 
ing Ansatz (|2^) for the staggered magnetization. A/ ~ 
L~f^l^ at i = , constitutes a further way of checking the 
2D Ising behaviour, since the critical exponents /3 = 1/8 
and z/ = 1 are involved. The data reported in Fig. ^ sistency for the 2D-Ising critical exponent ratios 



FIG. 6: Scaling of the staggered magnetization in the EA 
model with A,i = 0.01, for different t. The critical exponents 
/3 = 1/8 and u — 1 are those of the Ising universality class. 



give — 0.282(2). In the case of the weakest anisotropy 
Ap — 0.001, both the Binder's fourth cumulant and the 
staggered magnetization have not yet well entered the 
asymptotic scaling region for the lattice sizes considered, 
and ij cannot be reliably estimated by this technique. 

A further test of the universality class involves the lon- 
gitudinal staggered susceptibility x^^, Eq. (^): in this 
case the scaling Ansatz (GS) gives x 

(ij ~ L''/", with 

2D-Ising critical exponents 7 = 7/4 and ly — 1, as shown 
in Fig. 1^ for the case A^=0.001. The estimated criti- 
cal temperatures result <j(A^=0.01) = 0.2825(25) and 
tj(A^=0.001) = 0.2235(15), in full agreement with the 
above unbiased estimates. 

To summarize, in the case A^ = 0.01 we find con- 
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FIG. 7: Scaling of the longitudinal staggered susceptibility 
xT in the EA model with = 0.001, for different t. The 
critical exponents 7 = 7/4 and v = \ are those of the Ising 
universality class. 



FIG. 8: Scaling plot for the staggered magnetization in 
the EA model with A,^ = 0.01, for L = 16 (up triangles), 32 
(down triangles), 64 (diamonds), 128 (squares). The critical 
exponents /3 = 1/8 and v = \ are those of the Ising universal- 
ity class, and the critical temperature is taken as = 0.281 . 



and 7/1^, thus fully verifying the universality class. For 
= 0.001 the evidence, though limited to the matching 
of the estimates of ij obtained in Figs. ^ and |7[ is quite 
convincing. 

As a check that the magnetization and the staggered 
susceptibility have actually reached the asymptotic scal- 
ing regime with the considered lattice sizes, we have con- 
structed their scaling plots after Eq. (|2^), which are re- 
ported in Figs. H and ||. Data collapse for different lattice 
sizes is verified for the staggered susceptibility in the case 
A^ = 0.001 for L > 64, taking = 0.223, and a fortiori 
in the case A^ = 0.01; the staggered magnetization is in- 
stead seen to have entered the asymptotic scaling regime 
for L > 64 only in the case A^ = 0.01. 

From the above analysis a strong indication for the 
existence of an Ising phase transition is therefore given 
for both the considered anisotropics. Estimates of the 
critical temperature tj(A^) from the different criteria de- 
scribed in this section are summarized in Table |; all esti- 
mates are consistent, and amongst them we choose those 
realizing the best data collapse onto the universal scaling 
function in the scaling plots of staggered susceptibility 
and magnetization shown in Fig. |^ and ^j: the resulting 
values are t,(0.01) = 0.281(2) and (0.001) = 0.223(2). 
Such values will be indicated with a thin arrow in the 
following figures. 



B. Specific heat 

The specific heat of the Ising model is characterized 
by a sharp peak at the transition temperature. As the 
anisotropy decreases, a large bump, eventually coinciding 
with the bump of the isotropic model, grows on the right- 
hand side of the peak, which correspondingly moves to- 
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FIG. 9: Scaling plot for the longitudinal staggered suscepti- 
bility xr in the EA model with A^ = 0.001, for different L; 
symbols as in Fig. ^. The universal scaling function emerges 
from the overlap of the two solid lines. The critical exponents 
7 — 7/4 and v = 1 are those of the Ising universality class, 
and the critical temperature is taken as = 0.223 . 



TABLE I: 2D Ising transition temperature (A,, 
by FSS analysis and fit of critical behaviours. 



as obtained 



estimation method 


tj(O.Ol) 


(0.001) 


~ L 


0.2815(25) 


0.2225(15) 


Ui 0.6107 


0.280(3) 






0.282(2) 






0.2825(25) 


0.2235(15) 
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0.284(4) 
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FIG. 10: Specific heat of tlie EA modei vs t, for L = 64 (dia- 
monds), and 128 (squares); the dashed line represent the spe- 
cific heat of the isotropic model, as obtained by numerically 
deriving the internal energy QMC data of Ref. Arrows 
indicate the estimated critical temperature. 
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FIG. 11: Uniform susceptibility of the EA model for = 
0.01 and L — 64. Full diamonds: longitudinal; open dia- 
mondip-| transverse branch; stars: QMC data for the isotropic 
modeO. Solid and dashed lines are guides to the eye. The 
arrow indicates the estimated critical temperature. 




wards lower temperatures, meanwhile getting narrower. 
In Fig. |l^ we see that traces of an Ising-like peak emerg- 
ing from the isotropic curve can still be evidenced for 
both anisotropy values. Despite their being traces, we 
observe that they develop at the critical temperature as 
estimated above. These findings are in-good qualitative 
agreement with the experimental dataEil relative to the 
layered S = 1/2 antiferromagnet Cu(C5H5NO)6(BF4)2, 
which is supposed to have an anisotropy-driven transi- 
tion; similar behaviour is displayed by larger-spia. com- 
pounds whose anisotropy is known to be Ising-likel£3, such 
as K2NiF4 {S = 1) and K2MnF4 (S = 5/2). 



C. Uniform susceptibility 

At variance with the specific heat, where the 
anisotropic curves just slightly differ from the isotropic 
one, the uniform susceptibility undoubtedly shows an 
anisotropic behaviour: in Fig. ^l], where data relative 
to — 0.01 are shown, the transverse and longitudi- 
nal components, and x^^, separate from the isotropic 
curve att < 0.4, i.e., well above = 0.282. It is quite sur- 
prising that the Hamiltonian symmetry puts up so much 
resistance to the disordering effects of both quantum and 
thermal fluctuations: this means that the anisotropy, 
even as weak as those we are here considering, can never 
be neglected, and that there exists a temperature range, 
extending well above the transition (i.e., also out of the 
region where 2D correlations can trigger the onset of 3D 
long-range order), where genuinely 2D anisotropic be- 
haviour can be experimentally observed. 

The different temperature dependence of the trans- 
verse and longitudinal branches, with the former display- 



ing a minimum and the latter monotonically going to 
zero, is that expected for an EA antiferromagnet. This 
behaviour results from the anisotropy-induced spin or- 
dering, that makes the system more sensitive to the ap- 
plication of a transverse magnetic field, rather than of 
a longitudinal one. We observe that both the minimum 
of the in-plane component and the start of the rapid de- 
crease of the longitudinal one, are close to the transition: 
as such feature is peculiar to the Ising model, this re- 
sult gives further strength to the characterization of the 
transition as of Ising type. 

The two components of the uniform susceptibility are 
experimentally observable by means of conventional mag- 
netometry measurements: the above discussed deviations 
from the isotropic behaviour have been actually observeiji 
in severaLlayered compounds with S ^1: K2NiF4 



53 



Rb2NiF4£; 
BaMnF4E^ 



, BaNiF4Ei [S = 1), K2MnF4l^, Rb2MnF4S, 
(5 = 5/2). Such effects are here proved to be 
substantial also in S — 1/2 systems with a comparable 
anisotropy; unfortunately, to our knowledge, no clean ex- 
perimental realization of a 2D S* = 1/2 HAFM with small 
EA anisotropy is available yet. 



D. Staggered susceptibility 



The equal-time longitudinal and transverse staggered 
susceptibilities, and xf^, for — 0.01 are shown in 
Fig. [O- together with the susceptibility of the isotropic 
modeEj. Below the high-temperature region where the 
isotropic behaviour is reproduced, the two curves sepa- 
rate at i ~ 0.4, below which x^^ diverges more rapidly 
than in the isotropic case, while x^'^ stays finite and 
shows a maximum at about the transition temperature. 
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FIG. 12: Staggered susceptibility of the EA model for = 
0.01. Circles: longitudinal (bulk values); diamonds: trans- 
verse (L = 64); stars, lines and arrow as in Fig. |l^. 



The time-averaged susceptibilities display the same qual- 
itative behaviour, though their values are slightly differ- 
ent with respect to the equal-time case. 

As in the case of the uniform susceptibility, the ob- 
served behaviour is qualitatively suggestive of an Ising- 
like transition. Moreover, the analysis of longitudinal 
branch divergence gives a direct evidence of the Ising 
universality class, as well as an independent estimate of 
the critical temperature. For a 2D-Ising transition 
must display a power-law divergence, x^^ ~ |i — tj"''', 
with 7 = 7/4. In Fig. ^ we plot (x^/)'^^'' vs t for 
= 0.01, using data which are free of significant finite- 
size cor rectio ns, according to the criteria described in 
Section IV E . The power-law with the Ising exponent 
7 = 7/4 is evidently verified and the extrapolated criti- 
cal temperature is = 0.284(4), which agrees with the 
more accurate value obtained in Section IV A. As for 



the smaller anisotropy, — 0.001, the power-law di- 
vergence of x^^' could not be unambiguously detected for 
the considered lattice sizes. 



FIG. 13: Correlation length of the EA model for = 0.01. 
Symbols, lines and arrow as in Fig. |l^. 
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FIG. 14: Power-law critical behaviour of the longitudinal cor- 
relation length (diamonds) and of the longitudinal stag- 
gered susceptibility xT (open triangles) for = 0.01; solid 
and dashed line are linear fits of (^^^)~^'''' and {Xs^)~^^~' 
spectively. The critical exponents u = 1 and 7 = 7/4 are 
those of the Ising universality class. 



E. Correlation length 

Fig. |l^ shows the longitudinal and the transverse cor- 
relation lengths, and for A^ = 0.01. The two 
correlation lengths behave quite differently: the trans- 
verse branch, after having left the longitudinal one at a 
temperature t ~ 0.4, displays a maximum at the transi- 
tion, while the longitudinal branch diverges faster than 
in the isotropic model. Again, the overall behaviour is 
suggestive of a 2D Ising transition. 

The longitudinal antiferromagnetic correlation length 
is expected to display a power-law divergence ^ 
\t — t^\~'^, with V = I. One can capture this divergence 
by selecting a few points for at temperatures immedi- 



ately above t^, discarding those exceeding L/4, which are 
affected by finite-size saturation. This criterion is rein- 
forced by requiring the consistency of the estimates of 
obtained via the equal-time- and the time-averaged sus- 
ceptibilities: since both estimates converge to the same 
value in the thermodynamic limit, their agreement in- 
dicates that finite-size effects are under control. For 
A^ = 0.01, Fig. shows that is linear, with an 

extrapolated intercept i, = 0.283(6), in agreement with 
the value found via FSS analysis. 

The same observation is not possible for A^ = 0.001, 
as sizes larger than those here considered are required to 
approach the critical point of such model, while control- 
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FIG. 15: Longitudinal correlation length of the EA model 
with — 0.001, for L — 32 (down triangles), 64 (diamonds), 
128 (squares). Stars and lines as in Fig. hlL 



ling finite-size effects. 

We have also extracted the longitudinal correlation 
length in the vicinity of the critical point by fitting 
the equal-time correlator C^^(r), defined in Eqj-i|2l|), to 
a function due to Serena, Garcia and LevanyukO, 



F{x) = 



pl/2 



,1/4 



(31) 



properly symmetrized so as to take into account the pe- 
riodic boundary conditions, i.e., by 



C''{r) cx F{r/e') + F{{L - r)/S_'') 



(32) 



This function interpolates between the known asymptotic 
behaviours at r ^ and r — > oo of the Ising model. Well 
above the critical point we used the conventional fitting 
function for the isotropic antiferromagnetlfj 



F{x) 



(33) 



In the case = 0.001 good and stable fits are obtained 
even if the correlation length becomes comparable to (or 
even exceeds) the lattice size L: we can hence univocally 
define the fitted correlation length = ^|^. Moreover, 
the same kind of fitting procedure on the time-averaged 
correlator C^^(r), Eq. (pi]), gives consistent results. 

Notice that monotonically increases with L, 

bounded from above by the thermodynamic value; on the 
other hand, as is systematically smaller than , the 
latter is necessarily less sensitive to size finiteness. For 
this reason it is possible to observe in Fig. |l^ the clear 
deviation of from the isotropic model, due to its diver- 
gence at . To summarize, the sharp dependence of the 
longitudinal correlation length to small i^iiisotropies, |al- 
ready observed for S" = 5/2 in Rb2MnF4Ea and KFeF4Ea, 
is evidenced also for 5 = 1/2. 



V. EASY-PLANE MODEL AND BKT 
TRANSITION 

In this section we present our results relative to the 
EP model. We have used lattice sizes up to L = 200 and 
two anisotropy values: Aa ~ 0.02 (already considered 
in Ref. 10) and A^ = 0.001. These values are compa- 



rable with the experimentally estimated anisotropics of 
real compounds, amongst which several cupreous oxides 
such as La2Cu04, Sr2Cu02Cl2, aod Pr2Cu04, which are 
known to have an EP anisotropyS. 

The temperature range covered by our simulations is 
0.15|^S-.i ^ 0.90 : as suggested by previous calcula- 
tionslZlllJ this is the interval where we expect most of the 
peculiar features due to the anisotropy to be detectable. 
At higher temperatures the thermodynamic behaviour 
of the model does not differ from that of the isotropic 
one. On the other hand, finite-size limitations preclude 
the study of the very-low temperature region. To this 
respect, we recall that the correlation length of an EP 
model is expected to diverge exponentially as i — > tg^.^,; 
such fast divergence makes finite-size limitations more se- 
vere than in the EA case, where ^ diverges algebraically. 
On the whole, the BKT transition offers less robust evi- 
dences, both numerically and experimentally, due to its 
being a topological phase transition rather than a second- 
order one. 

In what follows we will refer to out-of-plane quanti- 
ties as those related to the hard z-axis, and to in-plane 
quantities as those related to the easy xy-plane. 



A. Finite-size scaling analysis 

The role of the staggered magnetization in the FSS 
analysis of the EA behaviour is somehow taken, i n the 
EP case, by the helicity modulus T, defined in Sec. Ill B. 



In the thermodynamic limit T is finite below and vanishes 
above the transition. When finite-size systems are con- 
sidered, the occurrence of a BKT transition is marked by 
the existence of a finite temperature below which T scales 
to a finite value, as suggested by Fig. ^ for Aa — 0.001. 

As for the value of the critical temperature, one knows 
that in the thermodynamic limit thcp^atio T/t gets the 
universal value 2/7r at the transitiontj. This behaviour 
is clearly detected in Fig. |l^, where the helicity modulus 
is shown vs temperature for different L: the slope of 
T{t) near the point where the line 2t/TT is crossed, gets 
larger for larger sizes, consistently with the occurrence of 
a jump in the thermodynamic limit. 

An upper bound to the BKT critical temperature can 
be hence given by looking at the temperature t where 
the scaling behaviour of T/t is most compatible with 
the expected asymptotic value 2/7r at criticality. From 
Fig. ^ we obtain t^^^{Ax=0.00l) < 0.180; the same 
procedure leads to t^^^{Ax=0.02) < 0.235. 

More accurate results are obtained by considering the 
critical scaling of T, Eq. (p9|). According to the procedure 
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FIG. 16: Scaling of T/f in the EP model for Aa = 0.001. The 
horizontal dashed line indicates the value 2/7r. 
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FIG. 17: Helicity modulus of the EP model for L = 16 (up 
triangles), 32 (down triangles), 64 (diamonds), 128 (squares), 
200 (circles). The dashed line is the function 2f/7r. 
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FIG. 18: Fitting parameter A vs t. The crossing point with 
the line A = \ gives an estimate of the critical temperature 
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FIG. 19: Scaling of t* with L for Aa = 0.02. The dashed 
line is a linear fit of the first three points, corresponding to 
L = 32, 64, 128. 



suggested in Ref. we assume the relation 
Ti(t) _ 2A{t) 1 



t 



1 



21og(L/Lo) 



(34) 



to hold in the vicinity of the transition; A{£) and io are 
then determined via a best-fit procedure and ^bkt identi- 
fied as the temperature where A{€) equals unity, as shown 
in Fig. |l^. The resulting estimates are t3^^(AA=0.02) = 
0.229(2) and ^3^,^ (Aa=0.001) = 0.172(5). In the case 
Aa = 0.001, this procedure is more uncertain: due to 
strong finite-size effects, T is seen to asymptotically scale 
just for L > 128 (to be compared with L > 32 in the 
case Aa = 0.02), so that the logarithmic fit can only be 
performed on two points {L = 128, 200) for each tem- 
perature. 

There is another way to exploit the data for the helicity 
modulus of a finite-size system, though we are not aware 
of such technique having been used by other authors be- 



fore. In Ref. |6^, Bramwell and Holdsworth found that in 
the classical 2D planar-rotator model on a finite size the 
ratio T/t gets the universal value 2/7r at a temperature 
t* > t^^^ , whose size dependence is given by 



t* 



4c(lnL)2 



(35) 



The above relation, determined by a renormalization- 
group based approach, is presented as a general prop- 
erty of BKT systems, though to our knowledge its va- 
lidity has never been checked for models others than the 
classical pure planar one. On the other hand, in the 
case Aa = 0.02 we can easily determine t* as a func- 
tion of L from Fig. |l^ and hence get Fig. ^ which 
shows that Eq. (|3^) holds even for weakly anisotropic, 
strongly quantum models; a linear fit of the scaling be- 
haviour of t* against (ln_L)~^ for L > 32 provides us 
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FIG. 20: Scaling of the in-plane staggered susceptibility Xs 
in the EP model with Aa = 0.001. 



TABLE II: BKT transition temperatures tgj^^(AA) as ob- 
tained by FSS analysis and fit of critical behaviour of 



estimation method 



,(0.02) 



W (0-001) 



asymptotic value of T 
A{t) = 1 
scaling of t* (L) 

~exp[ b^{t~t^^^y 



1/21 



< 0.235 
0.229(2) 
0.228(4) 
0.230(5) 
0.235(6) 



< 0.175 
0.172(5) 

0.180(5) 



with a rather accurate estimate of the critical tempera- 
ture, i (Aa=0.02) = 0.228(4). Moreover, the results 
of Ref. pll relate the coefficient c to the coefficient ap- 



pearing in Eq. ( |36|) in the form = n/^/c, and from 
the linear fit we obtain b^ — 0.96(9), in good agreement 
with the value obtained below by fitting the critical be- 
haviour of the correlation length. This remarkably shows 
that the predictions of Ref. ^ derived for the classical 
2D planar-rotator model, fully apply also to the quantum 
nearly-isotropic antiferromagnet we considered. 

Finally, an additional estimate of the BKT critical tem- 
perature is obtained by the in-plane staggered suscepti- 
bility, which is expected to scale at the transition as 
with rj = 1/4. Looking for the temperature where this 
scaling law is best verified, we obtain t^j^^ (Aa=0.02) = 
0.230(5) and t^^,^ (Aa=0.001) = 0.180(5), as shown in 
Fig. for the latter value. 

Although the identification of the BKT universality 
class is less complete than in the Ising case, substan- 
tial consistence between the different estimates of the 
critical temperature obtained by different predictions 
of the Kosterlitz-Thouless theory proves that the two 
anisotropic models display a BKT critical regime. The 
estimates for the critical temperature t ^^^ {A\) given 
in this section are summarized in Table O for the two 
anisotropics considered. Putting together these estimates 



FIG. 21: Specific heat of the EP model with Aa = 0.02 for 
1/ = 64 (diamonds) compared to QMC data fieji-feiie same 
modeltJ (triangles) and for the isotropic modelcJ'H (stars). 
The arrow indicates the estimated BKT temperature. Inset: 
zoom on the temperature region where a deviation is observed 
with respect to the isotropic case. 



we choose as reference values tg^^{A\—0.02) — 0.229(5) 
and t3^^(AA=0.001) = 0.175(10). Such values will be 
indicated by a thin arrow in the figures of the following 
sections. 



B. Specific heat 

The specific heat does not show large systematic devia- 
tions from the isotropic case within the resolution reached 
by the simulations for both anisotropics considered. Only 
a small temperature region, well above the estimated 
transition temperature, displays an anomaly in the form 
of a tiny peak, as shown in Fig. such peak is possibly 
reminiscent of the rounded peak shown by the specific 
heat of the.quantum S" = 1/2 XY model above its BKT 
transitionli^. We must however mentipj that, at variance 
with our results, previous QMC datallSl, also reported in 
Fig. ^ significantly deviate from the isotropic model. 
According to the generally low sensitivity shown by the 
specific heat to weak anisotropics, as seen for instance in 
the EA case, we find this result a bit unlikely. 



C. Uniform susceptibility 

As in the EA case, the uniform susceptibility reported 
in Fig. p2| shows strong evidences of the anisotropy. Mov- 
ing down from the high-temperature region, where the 
isotropic behaviour is reproduced, the in-plane ^ind 
out-of-plane uniform susceptibilities separate from 
each other and from the isotropic data. 

The in-plane component decreases faster than in the 
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FIG. 22: Uniform susceptibility of the EP model with Aa = 
0.02 and L = 64. Full diamonds: in-plane; open diampids: 
out-of-plane; stars: QMC data for the isotropic modeCJ'CJ. 
Solid and dashed lines are guides to the eye. The arrow indi- 
cates the estimated BKT temperature. 



isotropic case. At variance with the EA case, however, 
is not expected to vanish at t = 0, due to the contin- 
uous rotational symmetry of the ground state in the xy 
plane. Indeed, in a semiclassical picture, such symme- 
try allows the staggered magnetization to align along the 
in-plane axis perpendicular to the field, making possible 
the canted spin configuration with a finite ferromagnetic 
magnetization parallel to -and linear in- the field, so that 
X^^ stays finite. 

The out-of-plane susceptibility is instead enhanced 
with respect to the isotropic case, showing a minimum 
well above the transition. Such minimum marks the on- 
set of a completely different behaviour with respect to 
the isotropic model, entirely due to the presence of the 
small anisotropy. 



D. Staggered susceptibility and correlation length 

According to Kosterlitz theoryEl, in presence of a BKT 
transition the correlation length is expected to di- 
verge exponentially at finite temperature as 



= exp [b^it 



-1/2] 



(36) 



As for the = 0.02 model, this behaviour is in fact 
observed in Ref. |l^ where it is used to estimate the 
critical temperature. We use the estimates obtained 
by fitting the in-plane correlation function to Eqs. ( |3^ ) 
and (^3|). Discarding the values affected by finite-size 
saturation and thus considering only those satisfying 
^ we also observe the predicted behaviour: in 
particular, singling out the BKT critical region by succes- 
sively dropping points at high temperature until a stable 
fit is obtained, we obtain = 0.6(2), — 1.0(1), and 



FIG. 23: Correlation length of the EP model with Aa = 0.02. 
Circles: in-plane (bulk values); diamonds: out of plane (L = 
64); stars, solid line and arrow as in Fig. Ed. 
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FIG. 24: Staggered susceptibility of the EP model with Aa — 
0.02. Symbols as in Fig. |2^; stars, lines and arrow as in 



the estimate ig^T ~ 0.235(6), which agrees with the value 
found via FSS analysis. 

Furthermore, near criticality it is expected that the 
staggered in-plane susceptibility iwelated to the in-plane 
correlation length by the relationcS 



(37) 



where is a nonuniversal constant and ?/ = 1/4. 
By plotting ^^^^^ together with (xf^)^^^^~'''*i as done in 
Fig. one observes that this prediction also holds for 
the weakly-anisotropic quantum model; remarkably, the 
curve {x^^Y^'^'^~^'^ collapses onto the on a wide range 
of temperature so that if « 1: this property is not shared 
by, e.g., the classical planar rotator model or by the 2D 
quantum XY model. A closer look to the validity of the 
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FIG. 25: Critical behaviour of ^jf^ (full squares) compared to 

that of (xD^'''^""' with ri = 1/4 (open diamonds) for the EP 
model with Aa = 0.02. The solid Une is the BKT fit to the 
correlation length data. Inset: plot of InXa/ln^; the dashed 
line represents the expected BKT value 2 — rj = 1.75. 



scaling relation (^7|) can be obtained by plotting the ratio 



In Y^" 
In 



In if 



Inaj + fcai-iBKT)"'/' 



(38) 



which converges to the value 2 — 7] — 1.75 when t ^^kt' 
this is clearly shown by the data plotted in the inset of 

Fig. m. 

In the A\ = 0.001 case, neither ^f^ nor exhibit the 
expected BKT critical behaviour for the considered lat- 
tice sizes. The correlation length obtained by fitting the 
correlator C^^ to the function ( ^ ) is also not of much 
help. This suggests the in-plane correlation length to be- 
have as in the isotropic model up to relatively large values 
(^^xx ^ 100), and the same holds for the staggered sus- 
ceptibility. Such findings closely resemble those of neu- 
tron scattering experiments on very weakly anisotropic 
layered S—,= 1/2 compoiuiids, such as Sr2Cu02Cl2E3, 
La2Cu04E3 and Pr2Cu04Ej, that do not show signature 
of the existing anisotropy in the correlation length and 
static structure factor data. 

Both the out-of-plane staggered susceptibility and cor- 
relation length have a non-critical behaviour, with a max- 
imum well above the transition, at a temperature which 
roughly coincides with that of the minimum of the out- 
of-plane uniform susceptibility, marking the onset of an 
anisotropy-dominated regime. Both maxima are clearly 
decoupled from the transition temperature, at variance 
with the maximum of the transverse staggered suscep- 
tibility and correlation length observed in the EA case. 
To this respect we mention a definite disagreement with 
Rcf. |l^, where out-of-plane correlation length is conjec- 
tured to diverge exponentially when T — > 0. We show 
such conjecture to be wrong, as is clearly seen to 
saturate to a finite value. 



VI. PHASE DIAGRAM 

The detailed analysis presented above for the EA and 
EP models separately, is now composed to form the phase 
diagram t^ ^^^ vs A^_a in Fig. |2^, where our best esti- 
mates for the critical temperatures relative to the four 
models considered are shown, together with data from 
Refs. 1^ and [l^. Critical temperatures are seen to be 
strongly reduced with respect to the classical values, as 
given for instance in Ref. |7[ however, the diagram clearly 
suggests the critical temperatures to stay finite for any 
finite anisotropy, both in the EA and in the EP case, 
thus leading to the conclusion that quantum fluctuations 
cannot destroy the transition. 

We can actually see that the above conclusion is the 
consequence of a more general finding. If one numerically 
analyses ^^^^ (A^^x) finds that a logarithmic dependence 
is well consistent with our data, as shown in Fig. ^ Such 
dependencfipalready predicted by renormalization group 
techniquesBO, is rederived in Appendix |c| on the basis 
of a fully classical argument. It is found that 



ln(c/A^) 



and 



T 



ln(c'/AA) 



(39) 



(40) 



where c and c' are constants, while Pg is the spin stiff- 
ness of the classical isotropic model, entering the above 
expressions via the exponential divergence of its corre- 
lation length. The dominant effect of quantum fluctua- 
tions on such divergence is embodied in the spin stiffness 
renormalization; therefore, if the ordering process we are 
here observing is the same as in the classical case, we 
expect 47rps = 2.26 J, where the value Pg = 0.18 J has 
been taken for the renormalized = 1/2 isotropic spin 
stiffnes^. From the logarithmic fits of the quantum data 
we indeed find 2.22 and 2.49 as prefactors of the loga- 
rithm, which are remarkably near to the predicted value, 
despite the simplicity of the argument that led to it. 

For Aa = 0.02 and A^ — 0.01, where a direct compar- 
ison is possible, the critical temperatures are not fully 
consistent with the values given in Refs. ^ and |l^. We 
notice that the latter were estimated as free parameters 
of fitting functions for the critical behaviour of the sus- 
ceptibility and correlation length; the precision of this 
approach is hindered by the fact that the critical regime 
of both quantities was not always properly entered in the 
simulations of Refs. ^ and |l^, mainly due to technical 
limitations which are nowadays overcome. We therefore 
propose our data as more precise estimates. 



VII. CONCLUSIONS 

In this paper we have presented an extensive numer- 
ical study of thermodynamic and critical properties of 
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FIG. 26: Phase diagram of the S = 1/2 2d XXZ model on the 
square lattice for weak anisotropies. Full symbols are results 
of this work, open symbols are QMC data from Refs. |and0. 



is due to the fact that, in order to discriminate be- 
tween T — Q isotropic and finite-T anisotropic diver- 
gences one must come very close to the critical point of 
the anisotropic model, which is a non-trivial issue both 
numerically (due to severe finite-size effects) and exper- 
imentally (due to finite experimental resolution and in- 
tralayer coupling). 

As for the EP case, we underline that the considered 
values of anisotropy compare to that of several real com- 
pounds. On the other hand, we have clearly shown, for 
instance in Figs. |2^, and ^ that in the EP case 
traces of 2D anisotropic behaviour are detectable above 
the transition, due to the fact that some quantities dis- 
play either a minimum or a maximum in a temperature 
region well apart from tQ^Ti where the in-plane correla- 
tion length has not diverged yet, and experimental ob- 
servation should hence be more feasible. We therefore 
think that our results could constitute a sound basis for 
a possible experimental observation of genuinely 2D EP 
behaviour in real magnets. 



weakly anisotropic two-dimensional quantum antiferro- 
magnets described by the 2D S' = 1/2 XXZ model with 
both EA and EP anisotropy. Use has been made of the 
continuous-time QMC method based on the loop algo- 
rithm, implemented here for the first time also in the EA 
case, and on the worm algorithm, here reformulated as a 
variant of the loop algorithm. 

The general outcome of the numerical simulations is 
that the thermodynamics of 2D quantum antiferromag- 
nets is extremely sensitive to the presence of anisotropies 
of magnitude comparable to those of real compounds, i.e., 
as small as 10^"^ times the dominant isotropic coupling. 

In the models studied we see a finite temperature tran- 
sition to persist with clear signatures of Ising and BKT 
critical behaviour, in the EA and EP case, respectively; in 
the more anisotropic case (10~^) full consistency with the 
expected universality class is found. The most striking 
evidences of the presence of the exchange anisotropy are 
observed in the thermodynamic behaviour of correlation 
lengths and susceptibilities. Moreover, the dependence 
of the critical temperature on the anisotropy is found 
to be quantitatively consistent with the prediction rela- 
tive to the classical case, with properly renormalized pa- 
rameters. This tells us that quantum effects can neither 
destroy the transition, nor change the ordering mecha- 
nism responsible for the transition to occur and that our 
quantum models, despite having S* = 1/2 and very weak 
anisotropies, do actually behave as renormalized classical 
ones. Given the results of Ref. ^ for 5 > 1, we can say 
that this conclusion generally holds for quantum Heisen- 
berg antiferromagnets on the square lattice. 

As for the thermodynamic behaviour of the specific 
observablcs considered, we find all the non-diverging 
quantities to be highly sensitive to the anisotropy, while 
critical quantities show deviations, with respect to the 
isotropic case, which are generally harder to detect. This 
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APPENDIX A: ISING LIMIT OF THE WORM 
ALGORITHM 

We show that the Ising limit of the worm algorithm 
is nothing but the standard single-flip Metropolis algo- 
rithm, considering the simple case of an Ising antiferro- 
magnetic spin dimer. Independent transitions via single 
spin flip are only two: (i) |tT) lit) and {ii) its re- 
versal. Within the worm algorithm, the transition (i) 
occurs with probability 1, since the spin conflguration 
I IT) admits no bounce; this equals the Metropolis result, 
since a net energy gain J^/2 is obtained in the process 
(z). Starting from the configuration \[]) a bounce pro- 
cess is instead possible with a probability per unit time 
7r(3, h) — J'^/2; if the worm experiences even a single 
bounce, it will u-turn finding on its way back only pla- 
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quettes of type 1(2), on which it cannot bounce, until 
it reaches its tail leaving the configuration unchanged. 
Therefore the transition occurs with probability 1 —pih, 
where pib is the probability of experiencing at least one 
bounce. Since the number of bounces follows a Poisson 
distribution with parameter /37r(3, 5), we have that 



where H'^^^^ = 'H[^^'^ +rLf ' and J = (J^i, Ja)- Such 
expression stands as the direct quantum analogue of the 
classical expression as given in, e.g., Ref. ^ for the plane 
rotator model; in the limit of zero temperature it repro- 
duces the expression given by Ref. |65| . 

The QMC estimator @ can be obtained directly 
starting from (B5) in the case 5* = 1/2. In the 
continuous-time limit, the estimator for the bond ex- 
change operators f^^^ = [J^^ /2) SfSf_^^^ takes the 
form: 



■AXY) 



-/3J^/2 



fc>l 



fc! 



(Al) 



which is exactly the Metropolis transition probability for 
a process involving an energy increase /2. 



APPENDIX B: ESTIMATOR OF THE HELICITY 
MODULUS 

In this section we derive the QMC estimator (^7|) for 
the helicity modulus starting from its thermodynamic 
definition (p6|). The derivation is a finite-temperature 
generalization of the one given by Sandvik in Ref. in 
the context of Stochastic Series Expansion, to estimate 
the spin stiffness, i.e., at zero temperature. 

We start from the "twisted" XXZ hamiltonian, with 
the twist applied along the 1-direction of the lattice, as: 



(Bl) 



J^^ sin</.(5,-%^^ - S^,St+a.) + J^SfS. 



Z QZ QZ 



n2 



where di = (1,0) and Ti.2 is the term containing only 
bonds along the 2-direction, which remains unchanged. 

We expand the twisted Hamiltonian to second order in 
6 as: 



n{<p) = ni^ = 0) - ^ji 

where 



= -^r- X! {^i^i+d^ - ^i^i+di) (B3) 



is the 1-component of the spin current operator, and 



{XY) _ 



Carefully deriving the free energy with respect to the 
twist, i.e., taking care of the non-commutativity be- 
tween the 7t!((/) = 0), Ji and Ti^'^^^ operator, one ob- 
tains for the helicity modulus, averaged over the 1- and 
2— direction of the applied twist, the expression: 



T = - 



1 



2J^^L2 



{n^''''^) + jyr{JiO)-Jir))^ (B5) 



(B6) 



where iV^^^ is the number of (+)-kinks (of the type 
liiti+di) ITiii+di)) on the i,i + di bond, and anal- 
ogously for iVr^^. Therefore the estimator for the XY- 
energy takes the form: 



(H(^y)) = -bN+ + N-) 



(B7) 



where is the total number of (±)-kinks present 
in each co nfig uration. The current-current correlator 
present in (pa) can be decomposed into bond-pair con- 
tributions as follows: 



(B8) 

Taking into account the 5*= 1/2 constraint S^S^\a) = 0, 
one obtains that 

I' dr{ [t4 (0) - f r,^ (0)] [f + ^ (r) - f (r)] ) = 

^^4-^rdJ(^4-^7dJ 



1 



:^4 



(B9) 



0(0^) (B2) Putting together (§7|) and (|B9|) with (g5|) one obtains: 



^ = ^ {{Nt - N,-) + {N^ - N^-y 



(BIO) 



where (A^j^j) ^^'^ total number of ±-kinks on 1(2)- 
bonds. Since a (+)-kink and a (— )-kink affect a spin 
path crossing the kink by shifting it of a lattice spacing 
in opposite directions, we have that the spin-path wind- 
ing number can be expressed as 



^1(2) = (^1%) -^1(2))/^ 



(Bll) 



and this reduces the estimator ( |B1Q| ) to the expres- 
sion (|2l). 



APPENDIX C: CLASSICAL DESCRIPTION OF 
THE ORDERING MECHANISM 

We give here a sketchy description of the ordering 
mechanism in slightly anisotropic 2D magnets, referring 



20 



to the classical limit where the antiferro- and ferromag- 
netic cases are thermodynamically equivalent. 

EA case. We rewrite the Hamiltonian (|l|) in the clas- 
sical limit as 

i.d 

n' ^ ^ J2{s^s^^, + sysy^^ (ci) 

i,d 

where s = (cos 6* sin sin sin 0, cos 0) is a unitary clas- 
sical vector and Jd is the classical exchange constant. In 
the above form the Hamiltonian is written as the isotropic 
Heisenberg Hamiltonian plus, as long as <C 1, a 
small Ising perturbation 7i'. The isotropic Heisenberg 
model has-an exponentially divergent correlation length 
as T ^ OB, C « aTe^^'^s/^, where is the spin stiff- 
ness of the classical model. At very high temperatures, 
i.e., T ^ Jci, the spins are fully uncorrelated, so that 
the anisotropy has little effect. When T « J^i correla- 
tions set on and clusters of almost aligned spins form on 
the length scale Very roughly, one can imagine the 
spins of each cluster C to be fully aligned, so that the 
anisotropy term can be written as 

T-i' = '^"'i^^ cos(0i ~ (f>i+d) sm9i sm0i+d 

« JciA^e' ^sin^^c , (C2) 

c 

where 9c is the polar angle of the spin orientation on 
each cluster and border terms of order ^ are neglected. 
Hence, the anisotropy term creates an effective potential 
for the orientation of each correlated cluster that has two 
minima in 6'c' = and tt (up and down) separated by 



an energy barrier AE = JdA^^^. When ^ increases 
upon lowering T, the barrier becomes comparable to the 
thermal energy, so that the orientation of each cluster is 
confined to one side of the potential barrier: the system 
becomes Ising-like and a finite magnetization appears. 
Using the isotropic behaviour of ^, this happens when 

T^AE^ Jei A^ {aT e^'^Ps Z^) ^ . ((.3) 

solving this equation gives the critical temperature as in 
Eq. (g). 

The above simplified picture accounts for the loga- 
rithmic dependence of Tc upon the (small) anisotropy 
A^ that rwp,s earlier obtained via more sophisticated ap- 
proachesES; the logarithm appears to follow from the ex- 
ponential correlation length in the isotropic model. 

EP case. The above argument can be essentially 
rephrased, this time taking as perturbation of the 
isotropic Hamiltonian the term 

n'='^Y.'l<+d' (C4) 

which, in presence of clusters on the scale ^, becomes 

JciAa^' ^cos^^c • (C5) 
c 

The anisotropy potential has the minimum at 9c = Tr/2, 
i.e., for a cluster orientation on the xy-plane, and the 
well depth is AE — J^AxS.^. As in the EA case, the 
anisotropy becomes relevant once T is comparable to 
AE: the out-of-plane fluctuations are suppressed making 
the system effectively planar, so that vortex excitations 
appear and the BKT transition can take place. This 
roughly happens when Eq. (|4^) holds. 
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